Birds adapted to cold conditions show greater changes in range size related to past climatic oscillations than temperate birds

Investigation of ecological responses of species to past climate oscillations provides crucial information to understand the effects of global warming. In this work, we investigated how past climate changes affected the distribution of six bird species with different climatic requirements and migratory behaviours in the Western Palearctic and in Africa. Species Distribution Models and Marine Isotopic Stage (MIS) 2 fossil occurrences of selected species were employed to evaluate the relation between changes in range size and species climatic tolerances. The Last Glacial Maximum (LGM) range predictions, generally well supported by the MIS 2 fossil occurrences, suggest that cold-dwelling species considerably expanded their distribution in the LGM, experiencing more pronounced net changes in range size compared to temperate species. Overall, the thermal niche proves to be a key ecological trait for explaining the impact of climate change in species distributions. Thermal niche is linked to range size variations due to climatic oscillations, with cold-adapted species currently suffering a more striking range reduction compared to temperate species. This work also supports the persistence of Afro-Palearctic migrations during the LGM due to the presence of climatically suitable wintering areas in Africa even during glacial maxima.

Investigation of adaptive responses and distributional shifts of bird species to past climate changes provide crucial information to understand present and future effects of global warming and to adopt suitable conservation strategies. Past geographic distributions can be reconstructed with the help of the fossil record. GIS paleoclimatic layers and mathematical tools, such as Species Distribution Models (SDM), allow to project the current climatic requirements of species onto different past climatic scenarios, generating predictions of past distribution of species that are basically envelopes of climatic suitability, assuming niche conservatism over time [15][16][17][18][19] . Climatic envelopes are increasingly used to explore the potential distribution of species in the past and test evolutionary and biogeographical hypotheses. They are often integrated with molecular data to reconstruct detailed phylogeography and past population history of species 7,20-25. Many works explore LGM species population dynamics 23,[26][27][28][29][30][31] ; some of these are exclusively focused on bird taxa, aiming to clarify population dynamics in the refugia during climatic extremes, test niche stability, identify climate threats and effects to optimize current conservation efforts [20][21][22]24,[32][33][34][35][36][37][38][39][40][41][42] . Among these, very few use the information provided by the fossil record to assess palaeobiogeographic hypotheses or to calibrate (or validate) the LGM predictive models 20,22,39,41 and only two include the modelling of LGM wintering grounds of Afro-Palearctic migrants 20,39 .
The relationship between changes in range size (as effect of climate changes) and species thermal niche (climatic tolerance) has been theoretically investigated in mammalian species in the frame of past climatic oscillations, linking range shifts in Europe to niche optimum (cold-warm), and not to niche breadth 43 . In birds, LGM climatic ensembles. Pyrrhocorax graculus The LGM ensemble strongly overlaps with the MIS 2 fossil occurrences of P. graculus (Fig. 1a). In the present day, this sedentary species linked to cold climates survives with a relict distribution in the mountain areas of the mid-latitudes of the Palearctic, due to Holocene climate warming. In the LGM, as witnessed by models and MIS 2 fossil records, it experienced a considerable range expansion, spreading at lower altitudes due to the downward shifts of the mountain vegetation (Fig. 1a). The lowering of the upper limit of the tree-line favoured the expansion of the mountain pastures with rocky ravines and cliff faces that this species uses for feeding and nesting.
Bubo scandiacus The MIS 2 fossil occurrences are not well predicted by the LGM wintering distribution ensemble neither by the breeding distribution ensemble; in any case, the LGM wintering projection predicts the fossil distribution better than the breeding projection (Fig. 1b). The models predict for the LGM a clear southward shift of both the breeding and the wintering distribution of this species, currently restricted to the arctic and boreal areas. In detail, the LGM breeding range is completely shifted towards central and southern Europe with no overlapping with the current breeding areas, probably due to the expansion of the ice cap (Fig. 1b). According to the LGM models, this species was considerably more widespread than the present day, as approximately the whole Europe satisfied the conditions of climatic suitability of this species, if we consider both the breeding and the wintering projections.
Athene noctua The LGM predicted distribution of this species is well supported by the MIS 2 fossil occurrences, that are limited to Southern Europe (Italy, Spain, Portugal and Southern France) (Fig. 2a). A southward shift of the northern edges of the range due to the expansion of the ice cap is predicted by the LGM climatic envelopes, causing a reduction of the distribution area compared to the present distribution (Fig. 2a). The range reduction during LGM could have led to isolation of the different populations in the Mediterranean refugia (as during other glacial maxima 6 ) as this species is best adapted to temperate and warm climates and has a limited dispersal capacity 60,63,64 .
Perdix perdix The MIS 2 fossil occurrences of this species are not very well predicted by the LGM ensemble (Fig. 2b). The models, together with the fossil occurrences, indicate that this sedentary species linked to open environments could have undergone a southward shift of the northern edges of its current distribution during Table 1. Species thermal indexes of the six bird species object of the work. The STI is identified with the average temperature experienced by a species across its distribution (i.e., mean temperature) and defines the climatic tolerance of species. www.nature.com/scientificreports/  www.nature.com/scientificreports/ the LGM, causing a slight range reduction compared to the present one (Fig. 2b). This is possibly due to the expansion of the ice cap and the consequent presence of unsuitable climatic and environmental conditions at the northern latitudes. This species proves to be tolerant to cold climatic conditions as, even during MIS 2, it seems to persist up to latitude 50° N. Crex crex The LGM breeding ensemble is only partially supported by the MIS 2 fossil occurrences, that are very limited for this species (Fig. 3a). Based on the models, this species could have undergone a reduction of the northern edges of the breeding range due to the expansion of the ice cap, causing a range reduction compared to the present one (Fig. 3a). The Late Pleistocene fossil records of this species are very rare west to the Pyrenees 61,62 and the species is still absent from Western Europe, being mostly spread in Eastern and Central Europe. The models predict a similar distribution for the LGM. The LGM wintering African ensemble looks comparable to the current range, or even somewhat larger, but always limited to South Africa (Fig. 3a).
Coturnix coturnix The LGM breeding climatic envelopes support the MIS 2 fossil occurrences (Fig. 3b). The models suggest a reduction of the northern edges of the breeding ranges during the LGM compared to the present distribution, with a slight reduction of the distribution area. The LGM African wintering ensemble predicts a slight expansion of the conditions of climatic suitability in Africa during the LGM, compared to the present ones, except for the disappearing from the North African wintering areas (Fig. 3b). Nonetheless, we consider the predicted wintering presence of this species in the Sahara desert during the LGM as non-reliable, as during the LGM the Sahara was even larger than today in terms of latitudinal extent [65][66][67][68] . www.nature.com/scientificreports/  ). In the case of P. perdix, we assume that this could be linked to the anthropic impacts that this species suffered historically. The current populations, at least in Western Europe, are composed by a mixture of wild and farmed individuals that are released for conservation or hunting purposes. This aspect influenced boththe genetic diversity of the species, as the released individuals hybridize with the wild ones with a high risk of genetic introgressions, and the natural distribution of the species, due to the release of these commercial individuals in areas outside their "natural" climate envelope 69,70 . In addition, as this species nests on the ground in grasslands, it suffered massive declines during the XXth century due to habitat loss and degradation caused by agricultural intensification. Such an intense anthropic impact is not observed in the other investigated species. Therefore, we hypothesize that the current potential climatic niche and ecological requirements of P. perdix, modified by the anthropic impact and quantified using the present-day species distribution, might be different than the Pleistocene ones.
In the case of B. scandiacus, LGM map predictions do not overlap with the fossil occurrences especially in the case of the breeding areas. Here, our results indicate that the presence of this species in the European midlatitudes, with a concentration in Southern France, might be explained by wintering individuals, providing interesting hints about B. scandiacus phenology during the LGM (Fig. 1b). However, this prediction does not agree with 71 , where Southern France is included in the reconstructed LGM breeding range of this species. Furthermore, the species is reported to breed at the Middle Pleistocene fossil locality of l'Escal, in Southern France, based on the presence of sub-adult bones belonging to this species 61,72,73 . The wide distribution of this species during the LGM, evidenced by its fossil records, and more in general during the Late Pleistocene, contrasts with our narrow map predictions, and leads to hypothesize that the populations of this species had a wider climatic niche in the past. One of the possible reasons for this change, besides the climate warming, is the competition with the other Palearctic large owl, Bubo bubo. The distribution ranges of these two species widely overlapped during the Late Pleistocene, as witnessed by their fossil occurrences 61,62 . As suggested by some authors, B. scandiacus may have arisen in the Early Middle Pleistocene from Bubo ibericus sp. nov. in Southern Europe 74 . This would explain why B. scandiacus might not have always been tightly linked to glacial conditions. However, since the last late glacial, B. bubo might have increased its populations as a consequence of the climatic amelioration, forcing B. scandiacus to move to the high latitudes 71,74 . After the LGM B. scandiacus possibly changed also some morphological and ecological traits to adapt to the environmental constraints of inhabiting high latitudes. In the present day, white plumage colour might be due to the persistent snow cover in the Arctic and the habit of day hunting during summer might be related to the hours of daylight at high latitudes. However, as suggested by some authors, during the late Pleistocene, when the species was spread at the mid-latitudes, it might have been a crepuscular/night hunter throughout the whole year, and its plumage might have been darker (brownish) 71 .
The changes in the ecological requirements of both P. perdix and B. scandiacus could explain the scarce consistency between the LGM envelopes and the fossil occurrences in the Western Palearctic (Figs. 1b, 2b).

Range shifts.
Based on the STI of the six species here investigated, Pyrrhocorax graculus and Bubo scandiacus can be defined cold-dwelling species, whereas the others are more tolerant with higher average annual mean temperatures across their distribution (Table 1). During the LGM, the expansion of the ice caps represented a major constraint for the distribution of living organisms at high latitudes. Our results show a clear latitudinal shift for B. scandiacus (Fig. 1b) and a reduction of the northern edges of the range for Athene noctua, Perdix perdix, Crex crex and Coturnix coturnix (Figs. 2, 3). In P. graculus there seems to be no latitude shift, but it is just because this species never reached the higher latitudes during postglacial expansions, that were instead reached by B. scandiacus during warm phases as the current one. An altitudinal movement is instead clear in P. graculus, with shifts towards high altitudes during the warm phases and downwards during cold phases, following the alpine meadows where this species feeds (Fig. 1a). As a reference for the European and African topography, we provided the relevant elevation maps in Supplementary Figs. S4 and S5.
The models and the fossil occurrences together suggest that striking changes in range size between the present day and the LGM have occurred in the two cold-dwelling species B. scandiacus and P. graculus (Fig. 2). The former species expanded in the whole Europe during glacial times. It is reported as far south as Gibraltar (Gorham's Cave) in the Lateglacial 75 and in Southern Italy during the Lateglacial (Grotta Romanelli) 76 , MIS 2 (Grotta di Cardamone) 62,77 and MIS 3 (Ingarano) 78 . During the LGM, P. graculus followed the downward shift of the upper limit of the tree-line that, together with the retreat of the forest cover, led it to spread in a considerably larger area corresponding to approximately the whole non-glaciated Europe. The different ecological requirements of these two cold-adapted species led them to adopt different strategies to cope with warming climates such as during the Holocene: P. graculus moved towards higher elevations, whereas B. scandiacus moved towards higher latitudes, but both invariably reduced their range size, being considered "glacial relicts" 79 . These species were possibly also favored during glacial times by the expansion of open areas at the expense of forest cover, as they live in rocky areas (P. graculus) and open areas (B. scandiacus).
The models suggest, for the species adapted to temperate conditions (A. noctua, P. perdix, C. crex and C. coturnix), a slight reduction of the distribution during the LGM (compared to modern ranges), due to the southern retreat of the northern margins of their distribution owed to the expansion of the ice cap. The core of the mid-latitude distribution of these species remained unchanged throughout the last 20 ka (Figs. 2, 3). Despite a slight reduction in the size of their range due to major environmental constraints at the high latitudes, these www.nature.com/scientificreports/ species were possibly favored by the expansion of open areas typical of the cooler climatic phases, as they all exploit grasslands and/or steppe for their survival. Cold-dwelling species show large changes in range size between the LGM and the present day. Thus, we can infer that the thermal niche of the Eurasian birds might be a predictive trait of the potential range size variations due to climatic oscillations 43,44 . A visual examination of our map predictions show that changes in range size between the different climatic phases does not seem to be linked to the migratory behaviour, as birds with different migratory behaviours show similar variations in the range size whereas among the resident species we observed changes in range size of different magnitude. Likewise, the extent of range variations does not even seem to be related to the breadth of the thermal niche (see "Range" in Table 1). Evidence is here provided that cold adapted species are more threatened by the current climate warming and need more effective conservation measures, as their range is already considerably contracted and will suffer further reduction in the future. We expect the same difference, in reverse, for the species that are more adapted to warm climates (extreme range reductions during glacial periods and pronounced expansions during warm phases). These aspects should be investigated in further research.
As there is no standardized methodology to generate range map estimations, these are linked to the methodological choices of researchers. The LGM climatically suitable areas of P. graculus, C. crex and A. noctua modeled, respectively, by 36,39,41 , result smaller than those proposed in this work. All the three different approaches used algorithms that risk overfitting (e.g., GAMs) or Maxent with a hard threshold selection of 10th percentile of presences (forcing range maps to be small), resulting in maps that underestimate the real behaviour of these species. Our predictions have, in general, larger suitability values for Central and Northern Europe during the LGM. We believe that our maps are closer to reality, as these three species inhabit meadows and grasslands, which were widespread during the LGM. Furthermore, the northern distribution of these species during the LGM is documented by the fossil occurrences (see Figs. 1a, 2a, 3a and Supplementary Table S6).
LGM migrations. The LGM projections of the distribution of the two long-distance Afro-Palearctic migrants Coturnix coturnix and Crex crex indicate that climatically suitable areas existed both in Africa in the wintering grounds and in the Western Palearctic in the breeding grounds, providing further evidence for the persistence of Afro-Palearctic migration during the LGM (Fig. 3). These data agree with the view that the migratory behaviour is linked to an increase in climate seasonality and deterioration of winter conditions [80][81][82][83] . Also, our results suggest that there is no evidence of the loss of the migratory behaviour in Afro-Palearctic migrants during the LGM 39 and that, as reported in a recent research 84 , bird migration remained an important global phenomenon throughout the last 50,000 years. Our data therefore challenge the hypothesis that the migratory behaviour is mainly a phenomenon of interglacial periods (linked to the recolonization of northern de-glaciated areas in postglacial times) that considerably reduced during glacial phases [85][86][87][88] .
The paucity of the sub-Saharan African fossil record hinders the reconstruction of the presence of the Western Palearctic breeders in Africa in the past. To date, fossil evidence of the presence of C. coturnix and C. crex in Africa during the LGM is lacking. Nevertheless, the presence of Eurasian long-distance migrants in the Pleistocene fossil record of sub-Saharan Africa, together with the absence of medullary bone in these fossils, supports the existence of the Afro-Palearctic migrations during Pleistocene 39,89 . Among the LGM fossils of C. coturnix, two are located in Gorham's Cave (Gibraltar, Spain) and Ohalo II (Israel) (Fig. 3b), along two of the main Afro-Palearctic migration routes, possibly indicating that these routes were already used during the LGM.

Methods summary
Present distribution of the selected species, fossil occurrences, climatic data and species thermal indexes. The main ecological characteristics of the six species here investigated, together with their present-day distributions, are reported in the Supplementary Data S1. The data polygons representing the current distribution of the six species have been downloaded from the IUCN database [54][55][56][57][58][59] , with a single polygon for the sedentary species (Pyrrhocorax graculus, Athene noctua, Perdix perdix) and two polygons (one for the Palearctic breeding range and one for the wintering range) in the case of the migratory species 90 (Bubo scandiacus, Coturnix coturnix, Crex crex). Based on the IUCN data [54][55][56][57][58][59] , the parts of the ranges where the species are introduced have been excluded, whereas the areas where the migratory species are resident have been merged with both wintering and breeding ranges using R, version 4.0.3 (2020-10-10). The resulting polygons are what we used as "presence" polygons (present-day distributions). As the BRT (Boosted Regression Trees) models also require "absences" polygons, the latter have been created ad-hoc for each species. "Absence" polygons were built manually using QGIS, so that they closely surround the "presence" polygons (without using any specific distance to make the polygons). In this way, the current distribution limits of the investigated species can be predicted more accurately and the potential variables limiting species ranges can be identified. So, we aimed to select a calibrating data set that was able to codify, besides presences, also the real limiting factors for these species. "Presence" and "absence" polygons used for this paper can be downloaded from https:// www. dropb ox. com/ sh/ twij2 89ss2 gqmzk/ AABOY P4Tl6 9x4fp J4Zff 4A8xa? dl=0.
An extensive bibliographic search allowed us to identify a total of 648 Western Palearctic sites dated to the Late Pleistocene (from 130 to 11.7 ka BP) 91,92 that yielded fossil remains of at least one of the six species (identified with certainty to species level). We collected these fossil occurrences in a dataset. Multiple fossil occurrences of the same species from different stratigraphic units of the same deposit were counted separately due to their possible different age. Then, in order to identify those reliably dated to MIS 2 (29-14 ka BP), we thoroughly checked the age of each fossil occurrence, using the Radiocarbon Palaeolithic Europe database, v. 26 93 , which reports radiocarbon conventional ages BP, or the PACEA geo-referenced radiocarbon database 94  www.nature.com/scientificreports/ (95% CI) 96 . After this check, we found 48 fossil localities dating back to MIS 2, which are the ones that we used in this paper, and we georeferenced them. We report the list of the MIS 2 fossil occurrences of the six selected species in the Supplementary Table S6. It is worth mentioning that the relative abundance of fossil records in the western part of the Western Palearctic compared to the Eastern part is related to the higher number of fossil localities (mainly Palaeolithic sites) which have been investigated in the former area and shouldn't be regarded as an indication of the absence of the species from Eastern Europe. We also searched for African fossil records of the selected species dating back to MIS 2 in the literature, but none was found, due to the paucity of studies dealing with African fossil birds. Also, the Western Palearctic fossil occurrences of the migratory species do not necessarily indicate the breeding of those species in that specific locality, as they could also belong to individuals which were migrating. Only the finding of juvenile bones or medullary bone (calcium deposit in the bone cavity linked to egg-laying in female individuals) undoubtedly indicate the breeding of a species in a certain locality 97,98 . We downloaded rasters of climate data for the present-day and the LGM (21 ka BP) from ecoClimate 99,100 . In this dataset, eight different general circulation models (GCMs) were available (CCSM, CNRM, GISS, FGOALS, IPSL, MPI, MRI and MIROC, with a resolution of 0.5°) from the Couple Model Intercomparison Project (CMIP5) and Paleoclimate Modeling Intercomparison Project (PMIP3) working groups. We used all the eight simulations to model the present and LGM climatic envelopes of the six species.
To calculate the STI we coupled, using R (version 4.0.3), the IUCN polygons of the species' distribution with the CHELSA present-day climate raster and extracted the mean, median, minimum and maximum values of the Annual mean temperature experienced for each species across its distribution range, that is considered a valid estimate for the species' thermal niche 44 . For Afro-Palearctic migratory species (C. coturnix and C. crex), the STI has been calculated only for the Palearctic breeding ranges and using the Mean temperature of the warmest quarter instead that the Annual mean temperature.
The STIs of each species are provided in Table 1. The assessment of the climatic tolerance of each species has also been based on the general ecology of the species 60 . For instance, if we consider the two cold-dwelling species, P. graculus only lives in mountain areas above the tree-line whereas B. scandiacus lives at high latitudes mostly above the polar circle. On the other hand, all the other species are spread in Mediterranean areas also at the sea level, being adapted to temperate and even warm climatic conditions. These ecological traits are well supported by the Species Thermal Indexes here calculated (Table 1), that indicate P. graculus and B. scandiacus as the two species with the lowest values of average annual mean temperature experienced across their distribution.
Species distribution models. The present-day and the LGM climatic suitability of the six species have been modelled using BRT (Boosted Regression Trees). This is a flexible technique that allows to handle different types of predictor variables, express nonlinearities and interactions, and accommodate missing data. It consists in a boosting algorithm, that iteratively calls a regression-tree algorithm to construct a combination or ''ensemble'' of trees. Regression trees are good at selecting relevant variables and model interactions, whereas boosting combines large numbers of relatively simple tree models adaptively, overcoming the inaccuracies of single tree models and giving improved predictive performance [101][102][103][104] . All models were fitted in R, version 4.0.3 (2020-10-10), using packages Rpart (version 4.1-15) 105 and Caret (version 6.0-88) 106 .
The models were calibrated using all the 19 climatic variables available and the nowadays presence and absence polygons for each species, randomly sampling 300 species records from the presence polygons and 300 from the absence polygons. Absence polygons (generated in QGIS) were located surrounding the presence polygons, to avoid selectiong absence data far away from the observed distribution of the species. By doing this, we aim to detect the climatic variables that are actually limiting the current distribution range of the species. The selected data sets were successively randomly split in training data (80%) and testing data (20%). The selected model settings are lr (learning rate) of 0.1, tc (tree complexity) of 1 and bag fraction of 0.5. Evaluation of the model performance was calculated with the Area Under the Curve (AUC), ranging from 0.5 (random) to 1 (prefect prediction). All the values relative to evaluation and variable importance of each model are reported in the Supplementary Data S2. Models were then projected into the eight GCMs of the LGM cold climatic scenarios and the eight GCMs relative to the present-day climatic scenarios. The outputs are maps with each cell having an index of climatic suitability between 0 (no suitability) and 1 (suitability). Successively, in order to reduce the uncertainties related to individual model projections, we averaged the eight models related to the different GCMs in two ensemble forecasts for each species, one for the present-day and one for the LGM (Figs. 1, 2, 3). In the ensemble maps, the value of each cell in the grid indicates how many models (out of eight) predict the presence of the species in that given cell. The fossil occurrences were successively plotted in the LGM ensemble maps to test if the model outputs correctly predict the fossil distribution of the different species during the LGM. The values of the cells of the ensemble corresponding to each fossil occurrence have been reported in the Supplementary Data S3.

Data availability
All data generated or analysed during this study are included in this published article [and its supplementary information files]. www.nature.com/scientificreports/